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Abstract 

A probability density function for the variability of 
ensemble averaged spectral estimates from helicopter 
acoustic signals in Gaussian background noise was eval- 
uated. Numerical methods for calculating the density 
function and for determining confidence limits were ex- 
plored. Density functions were predicted for both syn- 
thesized and experimental data and compared with ob- 
served spectral estimate variability. 

Symbol List 


m 

envelope function 


even Fourier coefficient 

bm 

odd Fourier coefficient 

c(t) 

band-limited Gaussian composite 
function 

E[ } 

expectation operator 

/() 

chi-square limit function 

F 

window correction component 

»() 

noncentral chi-square limit function 

G 

window correction component 

H± 

window correction component 

Iu{) 

modified Bessel function of order v 

L 

lower confidence limit 

m 

Fourier index 

M 

number of Gauss-Legendre weights 
and abscissas 

n(t) 

band-limited Gaussian noise function 

N 

number of spectra in ensemble aver- 
age 

P 

probability density function 

P 

amplitude of sinusoidal signal 

9 

index for statistical moments 

Q 

window bias function 

r(t) 

band-limited Gaussian signal function 

R 

ratio of sinusoidal energy to Gaussian 
energy 


spectral density function 

•(*) 

periodic signal function 

t 

time 


T finite Fourier transform integration 

time 

Uj Gauss-Legendre abscissas 

U upper confidence limit 

V relative variability 

Wj Gauss-Legendre weights 

w(t) window function 

Wj window coefficients 

W confidence coefficient 

X a (t) even signal-plus-noise modulation 

function 

odd signal-plus-noise modulation 
function 

z(t) input function to energy detection 

system 

z a (t) even noise modulation function 

Zfc(t) odd noise modulation function 

y ensemble averaged spectral estimate 

y normalized ensemble averaged spec- 

tral estimate 

yi lower limit of integration 

y u upper limit of integration 

y normalized Bessel function argument 

z(t) output function from energy detection 

system 

a dummy parameter of degenerate 

hypergeometric function 

/ 3 area under tail of probability density 

function 

7 dummy parameter of degenerate 

hypergeometric function 

6i improvement to approximate root 

€ frequency mis-match index 

C dummy argument of degenerate 

hypergeometric function 

t jj dummy parameter of chi-square limit 

function 

6 dummy variable of integration 

i>(t) pure tone function 

fi mean 

/z normalized mean 

ft 2 2 nd central moment 


1 


ft2 normalized 2 nd central moment 

£ argument of chi-square limit function 


P 


T 


*0 


Xo 

x* 


0 


parameter for variety of functions and 
equations 

energy of random process or sinu- 
soidal signal 

approximation to inverse standard 
normal distribution 

time interval 

phase function 

degenerate hypergeometric function 

chi-square distribution 

noncentral chi-square distribution 

uniformly distributed phase random 
variable 

radial frequency 

time-average integration operator 


Subscripts: 


/ filtered function 

n narrow-band Gaussian noise 

r narrow-band Gaussian signal 

3 sinusoidal signal 

t random variable associated with 

function 

Introduction 

The variability of a single spectral estimate for a 
time series consisting of only band-limited Gaussian 
noise has been shown to follow a chi-square probability 
density function with two degrees of freedom [1-2]. It 
has also been shown that the variability is independent 
of the length of time over which data is collected. 
Increasing the temporal length of data improves the 
spectral resolution and reduces the bias of the estimate, 
but does not reduce the uncertainty inherent in a single 
spectral estimate. The method generally employed to 
reduce variability is to collect an ensemble of spectral 
estimates and then average the ensemble on an energy 
or mean power basis. The variability of the ensemble 
average is then given by a chi-square density function 


with 2N degrees of freedom where N is the number of 
independent spectra in the ensemble. 

When the time series for which a spectral estimate is 
desired substantially violates the assumption of Gaus- 
sian variability, using the chi-square density function to 
describe spectral uncertainty can be misleading. The 
acoustic signal generated" b^a helicopter in flight con- 
tains periodic impulsive noise that gives rise to harmon- 
ics of the rotor blade passage frequency in the power 
spectrum. Spectral estimates of helicopter noise at fre- 
quencies coinciding with blade passage harmonics can 
show significantly different variability than that pre- 
dicted by a chi-square function [3]. Quite simply, the 
periodic impulsive noise does not follow a Gaussian dis- 
tribution and, hence, a chi-square density is inadequate 
to describe the statistical behavior of its spectral esti- 
mate. 

The purpose of this paper is to show that the math- 
ematical methods necessary to describe the effect of pe- 
riodic impulsive components on the variability of spec- 
tral estimates already exist in the literature in books 
by Davenport and Root [4], Whalen [5], and Burdic 
[6] t as well as in a variety of scientific papers associated 
with signal detection methods [7-9]. Those methods are 
discussed very briefly in this paper and are presented 
in more detail in Appendix A. A Gaussian approxima- 
tion to the resultant statistical function is also provided. 
Numerical methods for evaluating the statistical func- 
tions, for integrating the functions, and for determining 
confidence limits are presented in Appendix B. A FOR- 
TRAN 77 computer program that evaluates theoretical 
confidence limits for chi-square distributions is listed in 
Appendix C. 

Finally, spectral estimate techniques will be dis- 
cussed briefly and the mathematical model will be com- 
pared with spectral estimates of real and synthesized 
data to examine its validity. Synthesized data will be 
used to most accurately represent the assumptions and 
approximations made in the development of the model 
and to permit an examination of the effects of spectral 
bias. Comparison with real data will include recordings 
of tiltrotor hovers at short range and helicopter long 
range approach flights. 

Mathematical Model 

Noncentral Chi-Square Distribution 

Power spectra may be obtained from time series data 
using any of three equivalent methods [2], One ap- 
proach is to perform a Fourier transform on the autocor- 
relation of the time series by the Blackman-Tukey pro- 
cedure. Another is to directly calculate a finite Fourier 
transform of the original time series data. A third way 
to obtain a power spectrum, and the approach to be 
modeled here, is to process the time series data with an 
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energy detection system [10] composed of three com- 
ponents that operate on the data in series as shown in 
Figure 1. The first component of the system is a narrow 
band-pass filter centered on the frequency of interest. 
The second component is a square-law detector that 
squares the output of the filter. The output of this 
detector is proportional to the instantaneous power in 
the filter band. The third component is an averager or 
integrator that is equivalent to a low-pass filter. The in- 
tegrator converts the instantaneous power output from 
the square-law detector into a time-averaged or mean 
power that is equivalent to a spectral estimate in the 
frequency band over the time of integration. 


Band-Pass 

Filter 


Square-Law 

Detector 


Integrator 


Figure 1. — Energy detector equivalent of spec- 
tral estimate. 


Correspondence of the energy detection system 
spectral method with the other spectral methods re- 
quires additional constraints. The first is that the width 
of the band-pass filter must be small compared to the 
center frequency of the filter. The second constraint 
is that the width of the low-pass filter must be equal 
to the width of the band-pass filter. If the integration 
time of the low-pass filter (which is the inverse of the 
filter width) is equal to the time length of data used 
with the other spectral methods, the spectral estimates 
are equivalent. As will be seen later, the integration 
time will be allowed to increase without limit for the 
purposes of the mathematical development. One con- 
sequence of this is the elimination of bias from the den- 
sity functions derived herein. Obviously, the practical 
calculation of spectral estimates requires a finite data 
set, but the mathematical ploy is intended to make the 
derivation tractable and will be tolerated for that pur- 
pose. The issue of bias will be raised again, briefly, in 
the discussion of spectral estimates in the data compar- 
ison section. 

When analyzing acoustic data from wind tunnel and 
outdoor experiments involving aircraft, it is convenient 
to distinguish between that part of the acoustic data 
arising from the operation of the aircraft and the back- 
ground, or ambient, part. The background will be re- 
ferred to as noise while that from the aircraft will be 
referred to as signal. The signal is further assumed to 
consist of a periodic and a random part. Certain ide- 
alizations of system behavior and acoustic signal and 
noise characteristics are made to simplify the analysis. 
The background noise is assumed to be a band-limited 
Gaussian process and the signal is assumed to be the 
sum of a periodic part and a band-limited Gaussian 
process that is independent of the noise process. The 


methods described by Davenport and Root [4] are uti- 
lized in the ensuing development. The input function 
to the energy detection system x(t) is given by the sum 
of the signal and noise as 

x(t) = s(t) + r(t) + n(t) 

where j(t) is the periodic part of the signal, r(t) is the 
Gaussian part of the signal and n(t) is the noise. 

The band-pass filter is assumed to pass all frequency 
components in the band without distortion and reject 
all other frequencies perfectly. The output from the 
filter ®/(t) is then given by 

*/(*) = */(<) + r /(0 + n /(0 

where J/(t) is a constant-amplitude sine wave with a 
constant frequency in the pass band, r^(t) is a narrow- 
band Gaussian random process associated with the sig- 
nal, and ny(t) is a narrow-band Gaussian random pro- 
cess associated with the background noise. Because the 
signal and noise random processes are independent, it 
is possible to define a composite narrow-band Gaussian 
random process 

c /W - r /(0 + n /(0 

where the variance of the composite process, tr*, is 
equal to the sum of the variances of the narrow-band 
processes associated with the signal, a 2 t and the back- 
ground noise, <r£. For consistency of notation the en- 
ergy associated with the sinusoidal part of the signal 
will be referred to as a 2 = P 2 f 2 where P is a constant 
amplitude. 

The integrator is assumed to be an ideal low-pass 
filter that passes all frequency components in the band 
without distortion and rejects all other frequencies per- 
fectly. The output from the integrator and, hence, from 
the energy detection system is given by z(t). The prob- 
ability density function of a sample of the output, 2 *, 
can be written (see Appendix A for derivation) 

*•>-(£) (^) 

which is referred to as a noncentral chi-square distribu- 
tion [5], This probability density function describes the 
variability inherent in a single spectral estimate when 
the time series contains a periodic part that is man- 
ifested as a harmonic in the frequency band. When 
there is no periodic part contained in the time series or 
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when the frequency band under consideration does not 
include a harmonic of the periodic signal, the density 
function is simply 



which is an exponential function or, equivalently, a 
chi-square function of two degrees of freedom. These 
density functions exactly correspond to results obtained 
by Burdic [6] for a matched filter detector in an active 
pulse-echo detection system. 

When an ensemble of N independent samples of the 
detection system output is averaged, the new random 
variable 


1 N 


l=l 


has a probability density function that can be ob- 
tained by analogy with the results for the multiple pulse 
matched filter detector given by Burdic [6] or by coor- 
dinate transformation of the normalized distribution of 
Whalen [5] 
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This probability density function describes the variabil- 
ity inherent in an ensemble averaged spectral estimate 
when the time series contains a periodic part that is 
manifested as a harmonic in the frequency band. When 
there is no periodic part contained in the time series or 
when the frequency band under consideration does not 
include a harmonic of the periodic signal, the density 
function is simply ^ 


p(y) = 


N 

alT(N) 



N - 1 

e ~Nv/°l 


which is a chi-square function of 2 N degrees of freedom. 
An interesting feature of this equation is that when 
there is no random component of acoustic signal in the 
time series, the c\ term can simply be replaced by a 2 . 
Consequently, the origins of the random components of 
the time series, at least insofar as they are Gaussian 
and independent, are unimportant to the shape of the 
density function. It is also clear that where random 
Gaussian components alone are present, their relative 
levels are unimportant and the absolute level has no 
effect on the relative variability about the mean. 


Gaussian Approximation 

When the number of spectra included in the average 
is large, the Central Limit Theorem can be invoked 
[11] to approximate the noncentral chi-square function 
with a Gaussian function. To do this, however, a mean 
and variance are required. The q th moment of the 
noncentral chi-square distribution is given by 

E lv q ] = / v q p(y)dy 
Jo 



where the degenerate hypergeometric function $ is 
given by [13] 


*(«,7i0 = 1 + (^) < + ^ (tJt) 



The mean of an ensemble average spectral estimate is 
therefore 


»=E[y] = * 2 c +a 2 . 

which indicates that, as mentioned in the introduction, 
allowing the integration length to grow without limit 
has allowed the derivation of a probability density func- 
tion that shows no effect of the bias associated with fi- 
nite data lengths. The variance of an ensemble average 
is given by 

A*2 = £[y 3 ] - M 2 = (^c + 2 < r c< T ?)/-N r 

from which a Gaussian probability density function can 
now be written 


that approximates the noncentral chi-square function 
for large sample Bize or when the ratio of periodic 
energy to Gaussian energy in the spectral band is great. 
Two parameters useful for quantifying the deviation of 
a statistical distribution from a Gaussian shape and, 
hence, the applicability of a Gaussian approximation 
are skewness 


7l = 


1 + 3 R 


V'tf [(l + 2fl) 3 / 2 


and (excess) Kurtosis 
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6 1 + 4 R 

72 “ * [(1 + 2 B ) 2 

where R — &g/&c ratio of periodic or tonal energy 

in the spectral band to the total Gaussian energy (signal 
and noise) in the band. 

By defining a variable that quantifies the relative 
variability of a spectral estimate about the expected 
value 


v = ^2 = I ^ + 2*2*? = I l + 2fl 

^ V^(^ + ff . 2 ) 2 \ N{1 + R) 2 

it is possible to show that when no periodic components 
are present in the time series, then 



which serves to reiterate the point made above. When 
there is a periodic component to the time series that is 
substantially greater than the total random component 
the relative variability is approximately 



and the relative variability decreases as the ratio of 
tonal energy to random energy increases. Clearly, the 
relative variability is always less when there is a tonal 
component present than when there is not. If the ratio 
of tonal energy to Gaussian energy increases without 
bound or if the number of samples in the ensemble av- 
erage increases without bound, both skewness and (ex- 
cess) Kurtosis approach zero and the Gaussian approx- 
imation approaches the noncentral chi-square density 
function 

Confidence Limits 

The confidence limits of a density function defined 
only for positive arguments are determined by integrat- 
ing under the tails 


and either 



P= f p(y)dy or 
Ju 


1-0= f p(y)dy 
Jo 


and solving for the lower confidence limit, L , and the 
upper confidence limit, V , where W = 1 — 2/3 is the 
confidence coefficient for a two-tailed interval. When 
there is no sinusoidal component in a spectral estimate, 
the density function is chi-square with 2 N degrees of 
freedom and it can be shown that the limits are given 
by L — £/N where ( is determined by solving 

i tf-i 

-(! + «-) = .-<£ L 

k=0 


and U = £/N where ( is determined by solving 


N-l 


1 * tk 

-o-HO^Efr 

fc =0 


The expected value of the normalized density function 
for noise only is just unity and the variance is given 
by l/N. When a sinusoidal component is present in 
a spectral estimate, the density function is noncentral 
chi-square with 2 N degrees of freedom and numerical 
integration must be performed. The details of the 
numerical integration scheme are given in Appendix 
B. Solving for the upper and lower confidence limits 
involves finding the root of each of the integrations. 
The numerical methods used for this process are also 
detailed in Appendix B. 

A comparison of the noncentral chi-square, X/ji 
and chi-square, Xo> distributions is shown in Figure 2 
below. For this example, both distributions have four 
degrees of freedom because they describe an average 
of two spectra. Both distributions have the same 
mean (or total energy) but the ratio of tonal energy 
to broadband energy is R = 10 (or 10 dB) for the 
noncentral distribution. The horizontal scale of the 
plot is normalized by the noise energy of the noncentral 
distribution so the total energy of both distributions is 
given by 1 + R. The noncentral distribution, denoted 
by the solid line, is noticably narrower than the central 
distribution, denoted by the dashed line. 

The upper and lower 80% confidence limits, ex- 
pressed in decibels, are asymmetrically spaced because 
of both the decibel scale and the asymmetry of the dis- 
tributions. The plot shows there is an 80% confidence 
that the actual spectral level is no more than 1.43 dB 
above and no more than 1.94 dB below the estimated 
level when the tone-to-noise ratio is 10 dB in the spec- 
tral band. This compares with an 80% confidence in- 
terval from 2.89 dB above to 5.75 dB below an estimate 
in a spectral band containing only broadband noise. As 
the tone-to-noise ratio approaches zero the shape of 
the noncentral chi-square distribution approaches the 
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shape of the chi-square distribution. As the tone-to- 
noise ratio increases the noncentral chi-square distri- 
bution becomes narrower and the confidence limits ap- 
proach the estimated spectral level. 



z 

Figure 2. — Comparison of noncentral chi- 
square, xfp and chi-square, Xo» distributions. 

A comparison of the noncentral chi-square distribu- 
tion, X/j> and a Gaussian approximation is shown in 
Figure 3 below. The noncentral chi-square distribution 
again has four degrees of freedom because it describes 
an average of two spectra, The tone-to-noise ratio is 
5 dB so R « 3.16. The Gaussian distribution has the 
same mean (or total energy) as the noncentral distribu- 
tion. The horizontal scale of the plot is again normal- 
ized by the noise energy of the noncentral distribution 
so the total energy of both distributions is again given 
by 1 -f R. The variance of the Gaussian distribution is 
given by (1 + 2 R)/N. The noncentral distribution, de- 
noted by the solid line, is noticably more skewed than 
the Gaussian distribution, denoted by the dashed line, 
which is symmetric. 



Figure 3. - Comparison of noncentral chi- 
square, X /{9 and Gaussian distributions. 


The plot shows there is 99.9% confidence that the 
actual spectral level is no more than 4.86 dB above and 
no more than 11.84 dB below the estimated level when 
the tone-to-noise ratio is 5 dB in the spectral band. 
The Gaussian distribution has an upper 99.9% confi- 
dence limit of 4 dB but an undefined lower confidence 
limit because the lower tail of the curve extends below 
zero. As either the tone-to-noise ratio or the number of 
spectra increases the Gaussian distribution more closely 
approximates the noncentral chi-square. 

A comparison of the 80% confidence limits from the 
noncentral chi-square distribution, and chi-square, 
Xqj distribution is shown in Figure 4 with a tone-to- 
noise ratio for the noncentral distribution of R ~ 10 
(or 10 dB). As the number of spectra included in a 
spectral estimate increases the confidence limits of both 
distributions tend to converge toward the estimate. The 
limits of the noncentral distribution, denoted by the 
solid line, are always within the limits of the central 
distribution, denoted by the dashed line. 



Figure 4. — Comparison of confidence limits from 
noncentral chi-square, X/p and chi-square, Xq, 
distributions. 


A comparison of the 99.9% confidence limits from 
the noncentral chi-square distribution, X/$» and Gaus- 
sian approximation is shown in Figure 5 with a tone- 
to-noise ratio of R « 3.16 (or 5 dB). As the number 
of spectra included in a spectral estimate increases the 
confidence limits of both distributions tend to converge 
toward the estimate and towards each other. The limits 
of the noncentral distribution, denoted by the solid line, 
are always above the limits of the central distribution, 
denoted by the dashed line. 

A comparison of the 90% confidence limits from the 
noncentral chi-square distribution, X/j> Gaussian 
approximation, the chi-square distribution, Xo> and its 
Gaussian approximation is shown in Figure 6 with a 
tone-to-noise ratio of R « 3.16 (or 5 dB). The Gaussian 
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approximations converge on their respective chi-square 
distributions as the number of spectra included in an 
average increases. The convergence tends to be more 
rapid for the noncentral chi-square distribution and its 
approximation. 



Figure 5. - Comparison of confidence limits from 
noncentral chi-square, and Gaussian ap- 

proximation to X/t* 
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Figure 6. - Comparison of confidence limits from 
various distributions. 


Graphs of 80% and 90% confidence limits are shown 
in Figures 12 and 13, respectively, for tone-to-noise 
ratios of -oo, 0, 5, 10, and 20 dB and for number of 
spectra ranging from 1 to 25. 


Comparison of Model with Data 


Spectral Estimates 

The development of the noncentral chi-square dis- 
tribution depends on the assumption that the inte- 
gration time of a Fourier transform is allowed to in- 
crease without limit. Practical evaluation of spectral 


estimates generally entails the Fast Fourier Transform 
(FFT) which is a discrete finite transform. Both the 
discretization process and the finite length data record 
introduce features to spectral estimates that are not 
considered in the development of the density function 
presented here. Since the model is to be compared with 
data that is necessarily finite, the effect of bias on a 
spectral estimate that is introduced by a finite trans- 
form must be considered. The bias in a spectral esti- 
mate of random noise has been shown to be roughly 
proportional to the second derivative of the spectrum 
and inversely proportional to the square of the period 
of integration [ 1 ]. For this derivation, though, only the 
bias introduced by a finite transform on a pure tone is 
considered. 

Let a pure tone of frequency u> 0 be represented by 


tf(t) = sin (u/ 0 t -f 

where ^ is a random phase distributed uniformly on 
[“ '*’>*’]• The magnitude of a spectral estimate at w e 
can then be written as 


QW = |©K )| 2 


where 


2 r T 

6(w e ) = — w(f)i?(f) coa(u/ e t)dt 

*2 [ T 

+ 7p J o w(f)t?(i) sin (u e t)dt 

is the finite Fourier transform of the windowed tone 
and w(t) is the window function. The expected value 
of Q(rp) is then given by 

Q = 2 “^ QWdi/) 

For a Dirichlet or rectangle window, w(t) = 1, with 
w e = u> 0 the expected value of the magnitude is Q = 1 
indicating no bias. The expected value of the magnitude 
of a spectral estimate for a windowed pure tone of a 
frequency different from the estimate frequency is then 
an estimate of bias introduced by the frequency mis- 
match and finite window length. 

For a discrete transform the spectral estimates are 
made only at discrete frequencies 


w e = kAw 

where Au> = 2 x/T is the resolution of an FFT, T is 
the integration time or length of a data record, and k 


7 




is a frequency index. The tone frequency may then be 
specified as 

w 0 = (J: 4- c) Au> 

where k 1 b the index that refers to the discrete frequency 
nearest the tone frequency and e = (u Q — u> c ) /Aw is 
an indicial measure of the difference between tone and 
estimate frequencies such that 


A number of common window functions [12] may be 
represented by 


The bias correction may now be written in terms of 
the window coefficients, Wj, estimate frequency index, 
and mis-match index, €, as 

Q = sine 2 (re)( — j 

with 

c= "»(2inb) + ‘ fl - 


w(t) — Wg — Wi COB (Au?t) 

4- W2 cos (2Awt) 

- W3 cos(3Awt) 

where various choices for the coefficients w^ define par- 
ticular windows. The coefficients may be chosen so that 
the window is power preserving insofar as the total en- 
ergy or mean power of a spectral estimate of a random 
process is unbiased on the average 



Coefficients of a number of window functions of this 
type [12] expressed in a power preserving form are 
shown in Table 1 below. 

Table 1* — Coefficients of selected windows. 


and 


B± 



2k + € \ 

(2 k + e) 2 -l) 

2k + € \ 

(2*+ f ) 2 -4; 

_2* + e_\ 
(2 k + e) 2 - 9 / 


For the case where the tone frequency and estimate 
frequency are the same, u e = u a or e = 0, a window 
introduces some bias with the correction given, for 
example, by 


Window 

w 0 

Wl 

w 2 

W 3 

Dirichlet 

1.0000 

0.0 

0.0 

0.0 

Hanning 

0.8165 

0.8165 

0.0 

0.0 

Hamming 

0.8566 

0.7297 

0.0 

0.0 

Blackman 

0.7610 

0.9060 

0.1450 

0.0 

Kaiser-Bessel 
4-sample, a = 3 

0.7463 

0.9236 

0.1823 

0.0226 

Blackman-Harris 
min 4-sample 

0.7063 

0.9614 

0.2782 

0.0230 


201ogi O [<2] = 201og lo [w 0 ] 

^ —1.76 dB for Hanning 
« —1.34 dB for Hamming 

In other words, the discrete finite transform spectral 
estimate of a pure tone with a frequency other than one 
of the discrete estimate frequencies will exhibit a peak 
lower then the magnitude of the tonal energy (bias) 
and will show the remaining energy from that tone 
spread across every discrete frequency of the estimate 
(leakage). The value Q gives the fraction of the original 
tonal energy that resides in the frequency bin u e nearest 
the tone frequency w Q . This fraction of the original 
tonal energy that remains in the spectral frequency bin 
is the value that should be used for the energy of the 
tonal signal, Qa\ rather than aj f in calculating the 
density function and evaluating confidence limits. 

The bias correction may be made for any arbitrary 
symmetric window in the same manner. The window 
correction component may be expressed as 
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(2 k + e) 2 -i*J 


where the window coefficients, Wj, are determined by 
discrete Fourier transform of the window. 

Synthesized Data 

Synthesized data were used to verify the mathemat- 
ical form of the noncentral chi-square density function 
and the numerical methods employed in its evaluation. 
Data records were constructed by adding white noise 
to sine waves. Each value of white noise for each data 
record was generated as an independent sample from a 
standard normal distribution by a commercial pseudo- 
random number generator and then scaled to the ap- 
propriate absolute level. The phase of the sine wave 
for each data record was generated as an independent 
sample from a uniform distribution on the interval (0,1) 
by a commercial pseudorandom number generator and 
then transformed to the interval (-t, t). Each data 
record was transformed by an FFT after a window was 
applied and the squared magnitude was then calculated 
to generate a spectral record. A spectral average was 
determined from the appropriate number of spectral 
records and the resulting spectral estimate at the fre- 
quency of interest was recorded for statistical analysis. 



Figure 7. - Comparison of biased and unbiased 
noncentral chi-Bquare distributions with an es- 
timate obtained from synthesized data. 


(a = 3) window was used, two-record spectral averages 
were calculated, and 1000 estimates at the frequency 
band nearest the tone, 19.457 Hz, were recorded. The 
histogram shows the approximate density function esti- 
mated from the synthesized data. The smooth curve 
shown by the solid line is the noncentral chi-square 
density function where the window bias correction was 
made. The smooth curve shown by the dashed line 
is the noncentral chi-square density function where no 
window bias correction was made. 

The bias corrected theoretical curve agrees very well 
with the histogram estimate from synthesized data. A 
parametric evaluation of the agreement between the- 
ory and synthesized data is summarized in Table 2 be- 
low. The parameters of interest, mean, standard devi- 
ation, skewness, and (excess) Kurtosis, are listed in the 
first column. The values for those parameters estimated 
from the distribution of spectral estimates are listed in 
the second column. The corresponding bias corrected 
theoretical parameters appear in the third column. The 
t statistic in the fourth column is for testing the null hy- 
pothesis that the estimated and theoretical parameters 
are equal. The critical value at the 50% level for 1000 
samples is t .59(999] = -6744 so the null hypothesis is not 
rejected for either the mean or the standard deviation. 
The critical value at the 50% level for infinite degrees of 
freedom is £.50(00] = -6742 so the null hypothesis is also 
not rejected for either the skewness or (excess) Kurtosis. 

Table 2. — Selected parameters of biased noncen- 
tral chi— square distribution. 


Parameter 

Estimate 

Theory 

t statistic 

Mean 

0.2182 

0.2193 

-.3791 

Std. Dev. 

0.0889 

0.0883 

0.2910 

Skewness 

0.6224 

0.6417 

-.2495 

Kurtosis 

0.4545 

0.5580 

-.6698 


Given the conflicting assumptions necessary to de- 
rive the model and correct for bias, as well as limitations 
imposed by the random number generator and the use 
of a discrete transform, it would seem that the non- 
central chi-square distribution appropriately describes 
the distribution of spectral estimates for a pure tone in 
Gaussian noise and that the numerical evaluation of the 
distribution is accurately accomplished. 

Tiltrotor Hover Data 


The results from one realization of this process are 
shown in Figure 7 above. A pure tone of 20 Hz at 
a level of 70 dB was added to noise at a level of 90 
dB. The sample rate was 2344 per second so a data 
record size of 2048 points gave a resolution of about 
1.14 Hz and yielded an average noise level of about 60 
dB in each frequency band. A 4-sample Kaiser-Bessel 


Acoustic data acquired during an XV-15 tiltrotor 
hover test were obtained to evaluate the applicability 
of the noncentral chi-square distribution to helicopter 
spectral estimates. Data from a single channel were 
converted to engineering units and spectra were calcu- 
lated from sequential segments of data with no overlap 
in the same manner as the synthesized data. Spectral 
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averages were determined from the appropriate number 
of spectral records and the resulting spectral estimate 
at the frequency of greatest magnitude was recorded for 
statistical analysis. 

Because the noise level and tone level were not 
known a priori, they were estimated from the spectra. 
The sample rate for these data was 24590 per second 
so a data record size of 8192 gave a resolution of about 
3 Hz. A 4-sample Kaiser-Bessel (a = 3) window was 
used, two-record spectral averages were calculated, and 
44 estimates were made from the limited amount of data 
available. An average spectrum of all the estimates is 
shown in Figure 8 below. The solid line is the average 
of all of the spectra while the dashed lines describe 
the envelope containing the 44 spectral averages of two 
spectra each. The peak spectral level occurs at 27.02 
Hz and the background noise per band was estimated 
from the spectral average curve to be about 56 dB at 
that frequency. The tone level in that frequency band 
was then found by subtracting the noise energy in the 
band from the total energy in the band. No correction 
was made for bias because the spectral levels used for 
this calculation were already biased. 
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Figure 8. — XV— 15 tiltrotor spectral average and 
envelope of spectra comprising the average. 


The estimated tone and noise levels were used to 
generate the noncentral chi-square density function for 
comparison with experimental data shown in Figure 9 
below. The histogram shows the approximate density 
function estimated from the hover data. The smooth 
curve shown by the solid line is the noncentral chi- 
square density function based on estimated tone and 
noise levels. The smooth curve shown by the dashed line 
is the chi-square density function for the total energy 
contained in the frequency band. 
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Figure 9. — Comparison of noncentral chi-square 
and chi-square distributions with an estimate 
obtained from XV-15 tiltrotor hover data. 

The noncentral chi-square curve is much narrower 
than the histogram estimate from the hover data. On 
the other hand, the chi-square curve is very much 
broader than the histogram. While the histogram 
seems somewhat closer in appearance to the noncentral 
theoretical distribution, the difference is striking when 
compared to the excellent agreement with synthesized 
data. 

An examination of assumptions used in deriving the 
theoretical distribution should indicate the source of 
disagreement between experiment and theory. Assump- 
tions made in deriving the theoretical function include 
requirements that the noise be Gaussian and that the 
tone, or the Fourier component of the periodic signal, 
maintain constant frequency and amplitude. The hover- 
ing tiltrotor was 100 feet above the ground at a horizon- 
tal distance of 500 feet from the microphone so propaga- 
tion induced variability is an unlikely culprit. The high 
tone-to— broadband ratio and the brevity of the test, 
30 seconds, make it unlikely that any non-Gaussian or 
nonstationary character of the noise would explain the 
disagreement. The blade passage frequency should be 
stable so it seems that source amplitude variations may 
account for the greater than expected variability. Ei- 
ther the source level or the distribution of energy among 
harmonics varied. 

Helicopter Approach Data 

Acoustic data acquired during a Sikorsky S-76 he- 
licopter approach were also obtained to evaluate the 
applicability of the noncentral chi-square distribution 
to helicopter spectral estimates. Data from a single 
channel were converted to engineering units and spec- 
tra were calculated from segments of data in the same 
manner as with the other data except for an overlap of 
about 35% between successive data segments. Spectral 
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overlapping was used to get as many independent esti- 
mates as possible in a short time interval. Harris [ 12 ] 
has suggested that transforms with this level of over- 
lap are essentially independent when good windows are 
used. Spectral averages were determined from the ap- 
propriate number of spectral records and the resulting 
spectral estimate at the frequency of the second blade 
passage harmonic was recorded for statistical analysis. 

Because the noise level and tone level were not 
known a priori, they were estimated from the spectra. 
The sample rate for these data was 2344 per second so a 
data record size of 2048 gave a resolution of about 1.14 
Hz. A 4-sample Kaiser-Bessel (a = 3) window was 
used, three-record spectral averages were calculated, 
and 35 estimates were made from a small fraction of 
the available data. The helicopter was approaching at 
a constant speed so the range was constantly shrinking 
and the tone level constantly increasing. The data 
were selected from a fairly short time interval at a 
relatively long range to minimize the relative effect 
of the changing range yet give sufficient data for a 
histogram estimate of the density function. 

An average spectrum of all the estimates from these 
limited data is shown in Figure 10 below. The solid 
line is the average of all of the spectra while the dashed 
lines describe the envelope containing the 35 spectral 
averages of three spectra each. The peak of the second 
harmonic of the blade passage frequency occurs at 48.07 
Hz and the background noise per band was estimated 
from the spectral average curve to be about 58 dB at 
that frequency. The tone level in that frequency band 
was then found by subtracting the noise energy in the 
band from the total energy in the band. No correction 
was made for bias because the spectral levels used for 
this calculation were already biased. 



Figure 10. — S— 78 helicopter spectral average and 
envelope of spectra comprising the average. 


The estimated tone and noise levels were used to 
generate the noncentral chi-square density function for 


comparison with experimental data shown in Figure 
11 below. The histogram shows the approximate den- 
sity function estimated from the approach data. The 
smooth curve shown by the solid line is the noncentral 
chi-square density function based on estimated tone 
and noise levels. The smooth curve shown by the dashed 
line is the chi-square density function for the total en- 
ergy contained in the frequency band. 

The noncentral chi-square curve is much narrower 
than the histogram estimate from the hover data. The 
chi-square curve shows much better agreement with the 
histogram and would appear to correctly describe the 
variability of the spectral estimates made from this data 
set. The helicopter approached the microphone from 
a considerable distance and it is virtually certain that 
propagation induced variability plays a significant role 
in the variability of the received signal. The frequency 
stability of the source is probably good but the highly 
directional nature of the sound from a helicopter rotor 
and variations in aircraft attitude during forward flight 
at low altitude may cause variations in the sound level 
emitted in the direction of the microphone. Examina- 
tion of the envelope of spectra in Figure 10 reveals that 
tone variability was approximately the same as back- 
ground variability. 
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Summary and Conclusions 

A probability density function for the variability 
of ensemble average spectral estimates from helicopter 
acoustic signals in Gaussian background noise was eval- 
uated. A brief summary of the noncentral chi-square 
and chi— square distributions and Gaussian approxima- 
tions to them was presented. The details of the develop- 
ment were presented in an appendix. Numerical meth- 
ods used to calculate the density function and determine 
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confidence limits were presented in another appendix. 
A FORTRAN program that implements the numerical 
methods to calculate confidence limits also appears in 
an appendix. 

Examples were plotted to show differences and sim- 
ilarities between the density functions. Plots were also 
presented to show differences between confidence lim- 
its versus the number of spectra included in an average 
for various confidence coefficients and tone-to-noise ra- 
tios. The noncentral chi-square density function was 
then compared with synthesized data that closely ap- 
proximated assumptions made in development of the 
model, with short range tiltrotor hover data, and with 
very long range helicopter approach data. 

The excellent agreement between synthesized data 
and theoretical curves indicates that the numerical 
methods and computer program worked as desired. The 
somewhat poorer agreement between the tiltrotor hover 
data and theoretical curves is likely an indication that 
the assumption of constant source level was violated. In 
this case the noncentral chi-square distribution was im- 
perfect but showed better agreement with the data than 
the chi-square distribution. The data from a helicopter 
approaching an observer showed very poor agreement 
with the noncentral chi-square distribution. The much 
better agreement of this data with the chi-square den- 
sity function is an indication that extremely variable 
tone levels, whether from source variability or propaga- 
tion effects, completely invalidate the use of the noncen- 
tral chi-square distribution. The noncentral chi-square 
distribution should give excellent agreement, however, 
over time scales where the tone levels do not vary sig- 
nificantly, as in wind tunnel measured acoustic data. 
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Appendix A: Mathematical Model 
Single Spectral Estimate 

This appendix contains the derivation of the probability density function of a sample of the output of an energy 
detection system with an input signal consisting of a pure tone ancTbroadband part along with broadband background 
noise. This derivation continues from the body of the paper. The sinusoidal part of the signal may be written as 

Sf(t) = PcoB(u,t 4- ip) 

where P is a constant amplitude, w, is the constant frequency of the component of the periodic signal contained 
within the pass band, and V is a random variable distributed uniformly on (0, 2ir). Over a time interval of length r, 
where 0 < t < r, the composite random process may be expressed as a Fourier series 

«/W = £ 

m=l 


/2xmt\ . / Y 

Om COS f — — 1 + bm SID ( — — 1 


where the coefficients 



«»(^)* 
l e ' (1),in (~ T*)* 


are Gaussian random variables that become uncorrelated as r increases without limit. The composite process may 
be rewritten as 


where 


Cf(t) = x a (t)cos(u,t) - Xb(t)sin(w § t) 


-(*) »n [(5=2 - «.) i] } 

*6(0 = X) [°™“n [(-=- i] - 6m cos 

so that the filter output may be written as the sum of the expressions for the sinusoidal part and the composite 
narrow-band process. Thus, 


Xf(t) = X a [t) cos[o/j$] - -Xfc(t)8in[ajjt] 


where 


X a (*) = Pcos(^) + x a (t) 


X b ($) = Psin(^) + 
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The random variables x a t and x^ t from x a (t) and Z(,(t), respectively, can be shown [4] to be independent and 
Gaussian with zero mean and a mean square of 


E 


c at 


= E 


B 6t] = f S{w)du =al=al + a* 
J J — OO 


where u is a dummy variable of integration, S(o>) is the spectral density, a\ is the variance of the composite narrow- 
band random process, a 2 is the variance of the narrow-band process associated with the signal, and a 2 is the 
variance of the narrow-band process associated with the background noise. For consistency of notation the energy 
associated with the sinusoidal part of the signal will be referred to as a 2 = P 2 / 2. The joint probability density 
function of the random variables x a t and x ^ is then 


from which the joint probability density function of the random variables X a t, and \j> can be determined: 

**-*■•« - h [(^5b) [(^) .-K— 

If the filter output is expressed in terms of a sinusoidal function ofw, with an envelope j4(t) and phase fat) that 
are slowly varying functions of time compared with cos(ci/,t): 


xj(t) = >l(t) cos[o/jt + fat)] 

then the envelope and phase may be written, respectively, as 

A(t)= [^ 0 2 (t) + X t a (t)] 1/J 

m [Si 

Since the only nonvanishing terms in the expressions for x a (t) and z&(t) are those that fall in the narrow filter band, 
the frequency components of the envelope and phase are confined to a similar band centered on zero frequency. The 
joint probability density function of the envelope and phase random variables can be written as 

p(At, 4>t, rp) = (^ i ) e-^ +p3 - 3AtPco <* t -^y 2 ^ 

The probability density function of the envelope can be determined by integrating over fa and (with 6 = fa — *0): 


p( A t) 

which can be written as 


-ft) 


e -(A’+P 2 )/ 2<r= 


f j_ r 2ir \±_ 

\ 2 ir J 0 2* U 


e {A t P COM B)t<r\ 


;J <iv| 


Mt) = (4|) .-W+'-V^Jb (^) 
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where I 0 () is the zero-order modified Bessel function. 

The output from the square-law detector is 

= 4 2 (i)cos 2 [w,t + ^(*)] = + i-A 2 (t)cos{2[w,i + ^(t)]} 

The integrator is assumed to be an ideal low-pass filter that passes all frequency components in the band without 
distortion and rejects all other frequencies perfectly. The low-pass filter width is the same as the band-pass filter 
width so that the square of the envelope function, with all its frequency components in the pass band, is passed 
by the integrator. The high-frequency component of the square-law detector output, with a frequency of 2u> 9 , is 
rejected because the band-pass filter width must be small compared with to satisfy the assumption that Cf{t) 
is a narrow— band Gaussian random process. The output from the integrator and, hence, from the energy detection 
system is given by z(t) where 


z(t) = (x 2 (*)) = l -A'{t) 

The probability density function of zt can be obtained by a simple transformation to give an expression 


P( z t) = 


(Jj) 



which is referred to as a noncentral chi-square distribution [5], This is the probability density function that describes 
the variability inherent in a single spectral estimate when the time series contains a periodic part that is manifested 
as a harmonic in the frequency band. 


Appendix B: Numerical Methods 
Density Function Evaluation 

Although the density functions appear to be relatively clear and concise, their numerical evaluation is difficult. 
When there is no periodic signal present the density function is the product of a constant, a power, and an exponential. 
In general the function exhibits a region where the power component dominates, a region where the exponential 
component dominates, and a region where the power and exponential components approximately balance. Because 
the magnitudes of the individual components can exceed the capacity of a computer to easily represent them while 
their product does not, the density function is evaluated numerically by expressing it in terms of an exponential 
only. Normalizing by the noise energy to simplify the expression gives 

p(y) _ e lnJV + (AT— l)ln(Ny)— JV S — Mr(W)) 

where y = y /ff 2 . The Gamma function is calculated by an asymptotic approximation [14] when N is sufficiently 
large. This expression is used even when a periodic signal is present if the ratio of tonal energy to broadband energy 
is less than 50 dB. 

When a periodic signal of significant magnitude is present the inclusion of the modified Bessel function enormously 
complicates the situation. Again normalizing by the noise energy to simplify the expression, the density function 
with a periodic signal is 


P(y) = w(!) ' e" Ar(t+,i) /iy- 1 (2JV v /«s) 

There are a variety of equivalent forma for a modified Bessel function of integer order but they each present difficulties 
for numerical evaluation. Four different forms, all from Abramowitz and Stegun [14], were chosen for evaluating the 
probability density function. The first equation uses a polynomial approximation to Iq that has a different form for 
each of two different regions 


p(y) * + 3.5156229/ + 3.0899424/ 

+ 1.2067492/ + 0.2659732/ 

+ 0.0360768/° + 0.0045813/ 2 ] 

for p < 1 and 


p(y) « [0.398942281 + 0.01328592p _1 + 0.00225319/T 2 

+ 0.00157565p~ 3 + 0.00916281/T 4 

- 0.02057706p _B + 0.02635537p“® 

- 0.01647633p -7 + 0.00392377p - ®]/ {/l Rjj 

for p > 1 where p = 8>/W/l 5. 

The second equation uses a polynomial approximation to 7 X that also has a different form for each of two different 
regions 


p(y) « 8ye _2 ( s+/i >[0.5 + 0.87890594p 2 + 0.51498869/ 

+ 0.15084934/ + 0.02658733/ 

+ 0.00301532/° + 0.0003241 1/ 2 ] 

for p < 1 and 
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p(y) * yW-R 3 c - 2 ('/»-' /S )’ [0.398942281 - 0.03988024P" 1 - 0.00362018p~ 2 

+ 0.00163801p -s - 0.01031555p -4 
+ 0.02282967p -5 - 0.028953 12p -6 
+ 0.01787654p -7 - 0.00420059p -8 ] 

for p > 1 where p = 16v^S/15. 

The third equation uses a uniform asymptotic expansion when both N > 2 and N -f log 10 (-R) > 3 

P(y) » e -iV (6+‘ R )+ ln ^+f(ln6-lnP)-tln(2iri/)-tln(l+y J )+^+l n ( l + ^® =j u *(p)/i/*) 

where i/ = JV — 1, y = 2(Nf w)y/Ry, p — 1/ ^/l y 2 , and the parameter tj takes the simple form 


1 = + y 2 +ln ^ if y < 2, 

1 + VI + y 3 

or the somewhat more complicated form 
and 

u fc+l(p) = ^P 2 (l - p 2 )w*(p) + | / (1 - 5p 2 )u fc (p)dp 

where uo(/>) = 1. 

The fourth, and final, equation uses an ascending series representation of the modified Bessel function such that 


K 

p(y) ~ Yl e^ +2&)lnAT+(7\r +Ac-l)ln y+fc In [r(fc+l)r(fc+AT)] 

*=0 


where if is determined when the K th term is sufficiently small compared with the sum. 

Integration 

A Gauss-Legendre numerical integration scheme takes the form 

where the M weights, wj, and abscissas, uj, are found on the interval (-1, 1) by approximating the roots of Legendre 
polynomials [15]. Because the probability density function sometimes exhibits significant values only over a finite 
interval, a further approximation can be made by replacing the extreme limits of integration so that the confidence 
limits are calculated by solving 

M 

Y^ w iP 

j=i 
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and 


1-/Jw 





for L and U . The lower limit of integration, X*, is arbitrarily set either to zero or to a value twelve standard 
deviations below the mean, whichever is greater 


L* = max (0, /i — 12^^) 

where the normalized mean and variance are given by /i = 1 -f R and = (1 + 2R)/N. 

Root Location 

Solving for the upper and lower confidence limits of a chi-square density function involves finding the root of 
each of two equations that exhibit the same mathematical form 


N-\ 


=«o = IT 


k= 0 


where tj — (1 4- W )( 2 for the lower limit and rj = (1 - W )/ 2 for the upper limit. Because both the first and second 
derivative of this function exist in analytic form 


and 



/"(■ o = / , (o(^ J -- 1 ) 

rapid convergence to the root can be achieved, given an acceptable initial guess, by using a refinement of Newton’s 
method [15] 


ti+1 


m 

fit*) 


A,- 


with 


A i = 1 + 


m rm 


where successive improvements are made to the approximate root until some convergence criterion is met. The 
function need only be evaluated once at each step and additional advantages are gained by observing that f'({) is 
the last term in the summation in f(() and also that only the ratio of f"(() to f'(() must be calculated. The second 
order correction term A< can cause instability so, in practice, its value is restricted to the range (0.1, 2.0). Because 
the summation is an approximation to the exponential function 


ff-l t k 

' £ »£|r 

k= 0 
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the function f(() is evaluated numerically by 


m = v - E 

k=o 

to minimize round-off difficulties. 

An initial guess for determining the root is provided by an approximation to the inverse chi-square function 
given by [14] 




where c is a rational approximation to the inverse standard normal distribution given by [14] 


? = ± 


/ 2.515517 + 0.802853p+ 0.010328p 2 \ 

V 1 + 1 432788p + 0.189269p 2 + 0.001308P 3 ) 


The sign of c depends on whether the tail is to the right or left and 


P = 



Because no simple expressions for the derivatives of the noncentral chi-square confidence limit functions exist, a 
secant method is used in which the first derivative term of Newton’s method is replaced by a secant approximation. 
The root of 


0 = 




which is the lower confidence limit, is found by making successive improvements to the approximate root 


L{+\ — Li + S t 


where 


= -$i_i 


9(Lj) 
g(Li) - 


until some convergence criterion is met. The correction term <5* can yield an approximate root outside the range for 
which the density function is defined so, in practice, its value is restricted to the range (- .99 L t {(X -- Li)/ 2). The 
same method is used for the upper confidence limit using the appropriate function with the exception that there are 
no restrictions to the range of 6+. 

An initial guess for the upper root is provided by the rational approximation to an inverse standard normal 
distribution, c, given above translated and scaled by mean and variance equal to those of the noncentral chi square 
distribution 


Ul = H + cv^2 
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An initial guess for the lower root is provided in nearly the same way when R is large, except that the value is not 
allowed to be too small 


L x = max(/i/10, n + cVmI) 

The approximation to the inverse chi-square function, given above, is used when R is less than 1/10. The first 
increment to the initial guess of the root, <5j, is arbitrarily set to /t 2 /10. 
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Appendix C: FORTRAN Program 

The FORTRAN program CHISQR calculates theoretical confidence limits for a chi-square distribution and a 
noncentral chi-square distribution. An example of program output is appended. 
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program chisqr 
C 

C 

C 

C CHISQR calculates confidence limits for chi-square and noncentral 
C chi-square probability density functions and compares them with 
C Gaussian approximations 
C 

C — 

C 

integer no 
parameter (no=51) 
real gla(no) ,glw (no) ,glx(no) 
integer io,n 

real aln , alo , arl , aru , csll , csul , del , eps , f n , f o , gcll , gcul 

real gncll,gncul ,nccsul,nccsll ,sdn,snn, snr ,snt 

real splr,spn,tt .vv.WjXp.xlmax.xlmin 

double precision dn,lnu,nu,snrl ,snrd,ldn,aO,al ,lngn 

common /gla/ gla 

common /glw/ glw 

common /glx/ glx 

common /dp/ n,snr .dnjlnUjnu.snrl.snrd.ldnjaOjal 

common /lngn/ lngn 

real xion.glt 

double precision lgamma 

eps=l . e-6 

sdn=12 . 

C 

C 

C 

C Calculate Gauss -Legendre integration weights and abscissas on (-1,1) 

C 

call gl(no,glx,glw) 
do io=l+(no+l)/2,no 

glx(io)=-glx(no+l-io) 
glw(io) -glw(no+l-io) 
enddo 
C 

C 

C 

C Enter the number of spectra, tone-to-noise ratio, and confidence level 
C 

write(6,699) 

write (6, 601) ’Enter the number of averages : ’ 

read(5,*) n 
C 

write (6, 600) 

write(6,601) ’Enter the tone-to-noise ratio (dB) : ’ 
read(5,*) splr 
C 
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write(6,600) 

write(6,601) 'Enter the confidence level (%) : 9 

read(5,*) w 
C 

C 

C 

C First guess from standard normal distribution inverse 
C 

aru= ,5*(1.+ .01*w) 
arl= .01*w) 
tt^sqrt (-2 . *alog(arl) ) 
xp=tt-(2 .515517+tt*( . 802853+tt* .010328) ) 
ft /(I . +tt*(l . 432788+tt* ( . 189269+tt* . 001308) ) ) 

C 

C 

C 

C Calculate tone and noise values 
C 

snr=10 . **( . l*splr) 
snrd=10 .dO** ( . Id0*splr) 
spn=l ,+snr 
C 

C Standard deviations of normal distributions 
C 

snn=l ./sqrt (float (n) ) 
snt=sqrt((l .+2. *snr)/n) 

C 

C Calculate preliminary values for EDSPDF 
C 

dn=dble(n) 
ldn=dlog(dn) 
a0=dn* (ldn-snrd) 
nu=dn-l .d0 

if (n . gt . 1 ) lnu=dlog (nu) 
snrl=dlog(snrd) 
al=2 .d0*ldn+snrl 
lngn=lgamma (n) 

C 

C Confidence limits for noise by chi-square distribution 
C 

csul=10.*alogl0(xion(arl t n) ) 
csll=10 . *alogl0(xion(aru,n) ) 

C 

C Confidence limits by Gaussian approximation 
C 

gncll=10 . *alogi 0 (amaxl ( 1 . 258925412e-10 , ( 1 . -xp*snt/spn) ) ) 
gncul=10 .*alogl0(l .+xp*snt/spn) 

gcll=10.*alogl0(amaxl(l . 258925412e-10 , (1 .-xp*snn))) 
gcul=10 . *alogi0( 1 . +xp*snn) 

C 

C — 
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c 

C Integration limits 
C 

xlmax=spn+sdn*snt 
xlmin=amaxl (0 . ,spn-sdn*snt) 

C 

c 

C 

C Search for confidence limits of noncentral chi-square using a 
C modified secant method for the lower limit and a secant method 
C for the upper limit 
C 

c 

C 

C First guess at lower limit 
C 

aln=amaxl( . l*spn , spn-xp*snt ) 

C 

C First guess if chi-square 
C 

if (splr .le .-10 . ) then 
vv=l . /(9 . ♦n) 

aln=amaxl (0 . , C ( 1 . -vv-xp*sqrt ( vv) ) **3) ) 
endif 
C 

C 

C 

C First function evaluation 
C 

f n=glt (xlmin , aln) -arl 
del= ( spn-aln) * . 05 
C 

C 

C 

C Repeat until there is an answer 
C 

1 continue 
alo=aln 
aln=alo+del 
C 

C Next function evaluation 
C 

f o=fn 

fn=glt (xlmin, aln) -arl 
C 

C Next delta 
C 

if (fo.eq.fn) goto 2 

del=amaxl (- . 99*aln , aminl ( . 5* (spn-aln) ,del*fn/(fo-f n) ) ) 
if (abs(del) .It .aln*eps) goto 2 
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goto 1 


C 

C Convergence 
C 

2 continue 

nccsll=10. *aloglO(aln/spn) 

C 

C 

C 

C First guess at upper limit 
C 

aln= spn+ xp * snt 
C 

C 

c 

C First function evaluation 
C 

fn=glt (xlmin , aln) -aru 
del=snt* . 1 
C 

C 

C 

C Repeat until there is an answer 
C 

3 continue 

alo=aln 

aln=alo+del 

C 

C Next function evaluation 
C 

f o=fn 

f n=glt (xlmin , aln) -aru 
C 

C Next delta 
C 

if (fo.eq.fn) goto 4 
del=del*fn/(f o-fn) 
if (abs(del) .It .aln*eps) goto 4 
goto 3 
C 

C Convergence 
C 

4 continue 

nccsul=10. *aloglO(aln/spn) 

C 

C 

C 

C Vrite results 
C 

write (6 ,602) ’Noise’ 
write(6,600) 
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write(6,603) * Chi-Square Gaussian 

write(6 , 603) 9 > f j 

write(6,604) 'Upper limit * ,csul , 9 dB',gcul,' dB' 
write(6 ,604) 'Lower limit csll , ' dB',gcll,' dB' 

C 

write(6,602) 'Tone + Boise, ',splr,' dB T/H Ratio' 
write(6 ,600) 

write(6,603) 'Boncentral Chi-Square Gaussian 

write(6 ,603) ' 

write(6,604) 'Upper Limit ' ,nccsul, * dB',gncul,' dB' 
write(6 ,604) 'Lower Limit ' ,nccsll , ' dB',gncll,' dB' 
write(6 , 699) 

C 

stop ' ' 

C 



c 

C Formats 
C 

600 format (x) 

601 format (x, a, $) 

602 format(////,x,a,2x,f5. l,a) 

603 f ormat (16x,a,a) 

604 formatU.a.ax.CSx.fO^.a.BxMSx.lO^.a)) 

699 format(/) 

C 

end 
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real function edspdf(z) 

C 



C 

C Foncentral chi-squared probability density function of order 2n 
C 



C 

real z 
integer n 
real snr 

double precision dn, lnu , nu , snr l t snrd f ldn f aO f al f lngn 
common /dp/ n,snr,dn,lnu,nu,snrl ,snrd,ldn,aO,ai 
common /lngn/ lngn 
integer k 

real qet a , qim , qin , qio , qip , qrho , qt , qt i , qt2 , qx , qzeta 

double precision eps 

double precision zz,zO 

double precision bO,bl 

double precision lnxi ,lnxim,pdf ,xi,lncf 

double precision s 

double precision zn,zs f eta, us ,den,z2 
double precision t ,tt , ttt ,tttt .ttttt 
eps=l .d-16 
C 



C 


C Evaluate PDF using polynomial approximation when N=1 


C 

C 


(Abramowitz and Stegun p 378) 


if (n.eq.l) then 
if (z.eq.O.) then 
edspdf = exp ( - snr ) 
else 

qeta=sqrt (snr) 
qzeta=sqrt(z) 
qx=2 . *qeta*qzeta 
qt=qx/3 .76 
if (qt.lt.l.) then 
qt2=qt*qt 
qio=l . 

+qt2*(3. 5156229 
+qt2* (3. 0899424 
+qt 2* (1.2067492 
+qt2* (0.2659732 
+qt2* (0.0360768 
+qt2*(0. 0045813)))))) 
edspdf =qio*exp(- ( snr+z) ) 
else 
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qti=l ,/qt 
qim=0. 39894228 

+qti*(0. 01328592 
+qti*(0. 00225319 
+qti*(-. 00157665 
+qti* (0.00916281 
+qti*(-. 02057706 
+qti*(0. 02635537 
+qti*(-. 01647633 
+qti*(0. 00392377)))))))) 
qrho=qzeta-qeta 

edspdf =qim*exp(-qrho*qrho)/sqrt (qx) 
endif 
endif 


return. 

endif 

C 



C 

C Evaluate PDF using polynomial approximation when N=2 

c (Abramowitz and Stegun p 378) 

C 

if (n.eq.2) then 
if (z.eq.0.) then 
edspdf =0. 
else 

qeta=sqrt (snr) 
qzeta=sqrt (z) 
qx=4 . *qeta*qzeta 
qt=qx/3 .75 
if (qt.lt.l.) then 


qt2=qt*qt 
qip=0 . 5 

+qt2*(0. 87890594 
+qt2*(0. 51498869 
+qt2* (0 . 15084934 
+qt2* (0 . 02658733 
+qt2*(0. 00301532 
+qt2*(0. 00032411)))))) 
edspdf =8 . *z*qip*exp(-2 . *(snr+z) ) 
else 


qti=l ,/qt 
qin=0. 39894228 

+qti*(- .03988024 
+qti*(- .00362018 
+qti*(0. 00163801 
+qti*(-. 01031555 
+qti*(0. 02282967 
+qti*(-. 02895312 
+qti*(0. 01787664 
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* +qti*(-. 00420059) ))))))) 

qrho=qzeta-qeta 

edspdf =sqrt (qzeta/ (snr*qeta))*qin*exp(-2,*qrho*qrho) 
endif 
end if 
return 
endif 
C 

C 

C 

C Evaluate PDF when N>2 
C 

if (z.eq.O.) then 
if (n.gt.l) then 
edspdf =0. 
else 

edspdf = exp ( - snr ) 
endif 
return 
endif 

zz=dble(z) 

if (snrd.gt . 1 ,d-5) then 

if (dn+dloglO(snrd) .gt . 3 .dO) then 
C 

C Uniform asymptotic approximation (Abramowitz and Stegun p 378) 

C 

zn=2 . dO* (n/nu) *dsqrt ( snrd*zz) 
z2=zn*zn 

if C zn . gt . 2 . dO ) then 

zs=zn*dsqrt (1 .dO+l .d0/z2) 
else 

zs=dmaxl(l ,dO,dsqrt(l .d0+z2) ) 
endif 
t=l .dO/zs 
tt=t*t 
ttt=t*tt 
tttt=tt*tt 
ttttt=tt*ttt 
if (zn.gt , 2 .dO) then 

eta=zs+dlog(l .dO/(l .dO/zn+dsqrt (1 .dO+1 .d0/z2))) 
else 

eta=zs+dlog(zn/(l .dO+dsqrtCl .d0+z2))) 
endif 
us=l .dO 
den=nu 

us=us+t*( . 125d0 

ft +tt * (- . 2083333333333333d0) ) /den 

den=den*nu 
us=us+tt* ( . 0703125d0 
* +tt*C-.4010416666666667d0 

ft +tt*( . 33420 138 8888 8 889d0) ) )/den 
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c 


ft 

ft 

ft 


ft 

ft 

ft 

ft 


ft 

ft 

ft 

ft 

ft 




den=den*nu 

us=us+ttt*( . 0732421875d0 
+tt*(-.8912109375d0 
+tt*(l .846462673611111d0 
+tt*(-1.02S812596450617d0))))/den 

den=den*nu 

us=us+tttt*( . 11215209960937Sd0 
+tt*(-2 . 3640869 140625d0 
+tt*(8 . 78912353515625d0 
+tt*(-ll .20700261622299d0 
+tt* (4 . 669584423426248d0) ) ) ) ) /den 
den=den*nu 

us=us+ttttt * ( . 2271080017089844d0 
+tt* (-7 . 368794359479632d0 
+tt* (42 . 53499874S38846d0 
+tt*(-91 .81824154324002d0 
+tt* ( 84 . 63621767460074d0 
+tt*(-28.21207255820025d0))))))/den 
den=den*nu 

us=us+ttt*ttt* ( . 572S0 142097473 14d0 
+tt* (-26 . 49143048695 155d0 
+tt* (218 . 19051 17442 116d0 
+tt * (-699 . 579627376 1326d0 
+tt* (1059 . 990452528d0 
+tt* (-765 . 252468 141 18 17d0 
+tt*(212.5701300392171d0)))))))/den 

xi=dn* ( (nu/ dn) *eta- (zz+snrd) ) 
xi=xi+ldn- . Sd0*dlog(zs)- . 5d0*lnu 
xi=xi+ . 5d0*nu*(dlog(zz)-snrl) 
xi=xi~0 ,918938533204673d0 
xi=xi+dlog(us) 
edspdf=sngl(dexp(xi) ) 
else 


C Ascending series (Abramowitz and Stegun p 375) 

C 

z0=dlog(zz) 

b0=a0+nu*z0-dn*zz 

bl=al+z0 

pdf=0.d0 

k=0 

lncf =lngn 
lnxim=bO-lncf 

1 continue 

lnxi=bO+bl*k-lncf 

if ((lnxi.lt .lnxim). and. (pdf .eq.O. d0)) goto 3 
if (lnxi.lt .-300. dO) goto 2 
xi=dexp(lnxi) 
pdf=pdf+xi 

if (xi.le.eps*pdf ) goto 3 

2 continue 


31 



k=k+l 

lncf =lncf+dlog(dble(k) ) +dlog(dble(n+k-l) ) 
goto 1 
3 continue 

edspdf =sngl (pdf ) 
endif 
else 
C 

C Broadband approximation when tone-to-broadband ratio is small 
C 

if (n.eq.l) then 

edspdf =sngl(dn*dexp(-zz) ) 
else 

if (z.eq.O.) then 
edspdf =0 . 
else 

s=zz*dn 

edspdf =sngl(dexp(ldn+nu*dlog(s)-(s+lngn) ) ) 
endif 
endif 
endif 
C 

return 

end 
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double precision function lgamma(n) 

C 

C 

c 

C Natural logarithm of the gamma function, G(n) 

C 

C 

C 

integer n 
integer i 

double precision nn,in,ins 
C 

c 

c 

C Evaluate gamma function (Abramowitz and Stegun p 257) 
C 

if (n. It. 100) then 
lgamma=0. 
if (n.gt.2) then 
do i=3,n 

lgamma=lgamma+dlog(i-i .d0) 
enddo 
endif 
else 

nn=dble(n) 
in=l . dO/nn 
ins=in*in 

lgamma= (nn- . 5d0) *dlog(nn)-nn 
lgamma=lgamma+0 . 918938533204673d0 
lgamma=lgamma+in* (+8 . 333333333333333d-2 
k +ins*(-2 . 777777777777778d-3 

k +ins* (+7 . 936507936507937d-4 

* +ins* (-5 . 952380952380952d-4) ) ) ) 

endif 
C 

return 

end 
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real function xion(a,n) 


C 

C 

c 

c 

c 

C Solve: 

C 

C 

c 

c 


k 

-x N-l x 

SUM for x given A and IT, return x/H. 

k=0 k! 


real a 
integer n 

double precision delt , eps ,f xio ,xin ,xio ,lfnk,d ,xiol 
real vv,t,xp 
integer k 
C 

C 

C 

C Modified second order Newton’s method when N>1 (Davis and Rabinowitz 
C p 114) (first estimate by Abramowitz and Stegun p 941 via p 933) 

C 

xin=-dlog(dble(a) ) 
if (n.gt.l) then 
vv=l./(9.*n) 
if (a. It. 0.5) then 

t=sqrt (-2 . *dlog(dble(a) ) ) 

xp=t- (2. 515517+t*(.802853+t*. 010328)) 

* /(I .+t*(l. 432788+t*(.189269+t*. 001308))) 
xin=n*(l .-vv+xp*sqrt(vv))**3 

else 

t=sqrt (-2 . *dlog( 1 .d0-a) ) 

xp=t- (2 .515517+t*( . 802853+t* . 010328) ) 

* /( 1 . +t * ( 1 . 432788+t+ ( . 189269+t* . 001308) ) ) 
xin=n*(l .-vv-xp*sqrt(vv))**3 

end if 
C 

eps=l .d-8 
1 continue 

xio=xin 
xiol=dlog(xio) 
fxio=dexp(-xio) 
lfnk=0.d0 
do k=l ,n-l 

lfnk=lfnk+dlog(dble(k) ) 
fxio=fxio+dexp(k*xiol-lfnk“Xio) 
enddo 

delt=(fxio-a)*dexp(xio+lfnk+(l-n)*xiol) 

d=dmaxl ( . ldO ,dminl (2 . dO , 1 . dO-delt * ( . SdO* ( (n-l )/xio-l . dO) ) ) ) 
xin=xio+d*delt 
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if (dabs (d*delt) .gt .dabs (xio*eps) ) goto 1 
endif 

xion=sngl (xin/n) 

return 

end 


35 



real function git (xmin,xmax) 


C 

C 

C 

C Adjust Gauss-Legendre abcissas and integrate the PDF 
C 

C 

C 

integer no 
parameter (no=Sl) 
real gla(no) ,glw (no) ,glx(no) 
common /gla/ gla 
common /glw/ glw 
common /glx/ glx 
real xmin,xmax 
integer io 
real hx,ta 
real edspdf 
C 

c 

C 

C Adjust Gauss-Legendre abcissas for interval (xmin.xmax) 
C 

hr= . 5* (xmax-xrain) 
ta=xmax+xmin 
do io=l , (no+l)/2 

gla( io)=xmin+hr*(l . +glx(io) ) 
enddo 

do io=l+(no+l)/2 ,no 

gla(io)=ta-gla(no+l-io) 

enddo 

C 

C Integrate the PDF over (xmin.xmax) 

C 

glt=0. 
do io=l,no 

glt=glt+glw(io) *edspdf (gla(io) ) 
enddo 

glt=hr*glt 

C 

return 

end 
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subroutine gl(n,x,w) 

C 



c 

C Calculate Gauss-Legendre integration abscissas and weights on (-1,1) 
C 



C 

real x(l),w(l) 
integer i,k,m,n 

real den , dp , dpn , d 1 , d2pn , d3pn , d4pn , e 1 
real fx,h,p,pk,pkml ,pkpl ,t ,tl ,u, v,xO 
C 



c 

C Find Gauss-Legendre integration abscissas and weights on (-1,1) 

C (Davis and Rabinowitz p 487) 

C 

m=(n+l)/2 
el=n*(n+l) 
do i=i,m 

t=(4+i-l)*3. 141 5926536/ (4*n+2) 
xO=(l.-(l.-l. /n)/(8 ,*n*n))*cos(t) 
pkml=l . 
pk=xO 
do k=2,n 
t l=xO*pk 

pkpl=tl-pkml-(tl-pkml)/k+tl 

pkml=pk 

pk=pkpl 

enddo 

den=l .-xO*xO 
dl=n* (pkml-xO*pk) 
dpn=dl/den 

d2pn=(2 . *xO*dpn-el*pk)/den 

d3pn= (4 . *x0*d2pn+(2 . -el ) *dpn) /den 

d4pn= ( 6 . *x0*d3pn+ ( 6 . - e 1 ) *d2pn ) /den 

u=pk/dpn 

v=d2pn/dpn 

h=-u* ( 1 . + . 5*u*(v+u* (v* v-d3pn/ (3 . *dpn) ) ) ) 
p=pk+h*(dpn+ . 5*h* (d2pn+h/3 . * (d3pn+ . 25*h*d4pn) ) ) 
dp=dpn+h* (d2pn+ . 5*h*(d3pn+h*d4pn/3 . ) ) 
h=h-p/dp 
x(i)=xO+h 

f x=dl-h*el*(pk+ . 5*h*(dpn+h/3 . *(d2pn+ . 25*h* (d3pn+ .2*h*d4pn) ) ) ) 
w(i)=2 . *(1 ,-x(i)*x(i) )/(lx*fx) 
enddo 

if (m+m.gt.n) x(m)=0. 

return 

end 
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Example run of FORTRAN program CHISQR: 


Enter the number of averages : 5 

Enter the tone-to-noise ratio (dB) : 5 
Enter the confidence level (%) : 80 


Noise 



Chi-Square 

Gaussian 

Upper limit 

2.0377 dB 

1.9679 dB 

Lower limit 

-3.1290 dB 

-3.6978 dB 


Tone + Noise, 5.0 dB T/N Ratio 


Noncentral Chi-Square 


Upper Limit 1.4144 dB 

Lower Limit -1.9070 dB 


Gaussian 


1.3758 dB 
-2.0253 dB 


STOP 
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